Vortices in a rotating BEC under extreme elongation 
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We investigate a non-axisymmetric rotating BEC in the limit of rotation frequency for which the 
BEC transforms into a quasi-one-dimensional system. We compute the vortex lattice wavefunction 
by minimizing the Gross-Pitaevskii energy functional in the lowest Landau level approximation for 
different confinement potentials. The condensate typically presents a changing number of vortex 
rows as a function of the interaction strength or rotation-confinement ratio. More specifically, the 
vortex lattices can be classified into two classes according to their symmetry with respect to the 
longitudinal axis. These two classes correspond to different local minima of the energy functional 
and evolve independently as a function of the various parameters. 
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INTRODUCTION 

In recent years several groups have studied theo- 
retically and experimentally vortices in rotating Bose- 
Einstein condensates (BEC's). Butts and Rokhsarp] 
were the first in, with the help of the Gross-Pitaevskii 
(GP) functional, calculating the wave function of a ro- 
tating Bose-Einstein condensate using a lowest Landau 
level (LLL) approximation^. Interestingly, a very sim- 
ilar approach was being developed at the same time 
to study a class of closely related systems: Vortices in 
superconducting disks[2j. The theoretical idea for in- 
ducing the mechanical rotation and coupling of internal 
states using an electromagnetic field (phase imprinting) 
in the condensate was proposed later by Williams and 
Holland|| and the experimental confirmation was given 
shortly afterwards by the JIL A group 5] . Arrays of few 
singly quantized vortices were produced later with opti- 
cal spoon stirring method0|. After these works, it was 
established that in a rotating BEC the energetically fa- 
vorable state was a triangular lattice of vortices similar 
to the Abrikosov lattice in superconductors Ho, based 
on this, showed that a three dimensional BEC with a big 
number of vortices can be described like a quantum-Hall 
system in the LLL0. However, in the "ultrafast" rotat- 
ing limit, where the rotation frequency Q is close to the 
trapping frequency w, the GP mean field approximation 
fails 9]. The study of vortices in the ultrafast limit can 
still be done including a quartic radial potential in the GP 
functional which further confines the bosons and allows a 
mean-field description of the BEC above the critical ro- 
tation frequency |l(J. This kind of trapping potential has 
allowed to find experimental evidences of the melting of 
the vortex lattice predicted in previous works 



THE GROSS-PITAEVSKII FUNCTIONAL IN 
THE LANDAU GAUGE 



The vortex lattice state of a BEC in symmetric po- 
tential traps has been well studied with Gross-Pitaevskii 



and Thomas- Fermi formalisms [l2j. Recently, however, 
the behavior of the vortex lattice in an anisotropic poten- 
tial trap has caught the attention of investigators|13| and 
interest to look for unsual vortex structures in BEC's un- 
der different conditions has risen Anisotropic BEC's 
are modeled adding an asymmetric potential to the usual 
harmonic confining potential which elongates the BEC. 
To find the wave function of a rotating and quasi-two- 
dimensional BEC|l5l| confined by a symmetric harmonic 
potential one minimizes the GP functionally 



E[* ] = J dV (E k + E h + E t + E r ) 



(1) 



where E k 

mu> a (x +y ) 



= £-|V* | 

2m 1 u 1 

l^ol 2 is the harmonic term, 



is 



is the kinetic term, Eh 
= §l*o| 4 

the interaction term, and E r = — ^^fl ■ L^o is the rotat- 
ing term, with f2 the rotating angular frequency, and L 
the angular momentum operator. We assume through- 
out an oblate BEC of witdh a z in the z direction which 
keeps m the ground state for this direction. We take 
il = ilz so the rotating term becomes 



d_ _ d_ 

dy ^ dx 



(2) 



We now consider the anisotropic trap which produces the 
elongation of the BEC. This trap is modeled by adding 



dV 



mA 2 (x 2 



-A*o| 2 , 



(3) 



with A varying between and loq 0,0]. The total trap- 
ping potential is thus V to t = i j[x 2 (uo+A 2 )+i/ 2 (ljq-A 2 )]. 
When A grows, the system becomes more confined in 
the x direction and more elongated in the y axis. These 
kind of systems have been studied theoretically in pre- 
vious works for slightly elongated rotating traps [lil HI| 
and, recently, for extreme elongation as well|20|. Setting 
A 2 = ujg — il 2 , the system losses the y confinement and 
the BEC resembles a channel of length I, which tends to 
infinity, with parabolic confinement in the x direction. 
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Using the gauge transformation "to = r" 102 ' 1 '!^ the 
energy functional can be written as 



\ 2m 2 ay 

(4) 

where = ma; 2 (f2 2 + w 2 ) is the effective trapping 

potential. The relation between A and Q defines the ex- 
istence of only one rotating critical frequency for a given 
anisotropy for which the BEC is just deconfined in y. 

We use the oscillator length ao = yj and the har- 
monic ground state energy E = as length and en- 
ergy units, respectively. We propose an expansion of the 
wavefunction f as a linear combination of lowest level 
eigenfunctions, <f>k, of the kinetic energy operator, the 
rotating term, and the trapping potential: 



(5) 




FIG. 1: Energy as a function of confinement for the two fam- 
ilies of configurations: even (dashed-dotted) and odd (solid). 
The ground state corresponds to the lower envelope. 



By symmetry we propose eigenfunctions of the form 

4> k =Ae* k y Xk (x), (6) 

where A is a normalization constant. The one- 
dimensional differential equation for Xk becomes 



dx 2 



Xk 



>2Rx 



2 n 

-= — / 



Xk = 0, (7) 



with e* k = e k 



fc 2 



i+ 



, efe being the eigenenergies of 



the eigenfunctions <pk, and R = 1 + We find the set 
of normalized functions Xk numerically and analytically 

and set the normalization constant A to — != with V — 

vv 

aola z . In the rotating frame the functions 4>k are similar 
to the solutions of the Schrodinger equation for a charged 
particle in presence of a magnetic field described by a 
vector potential in the Landau gauge 21J. An analysis of 
the analytical solution of the differential equation can be 
found in the Appendix I. 

Now the GP functional becomes 



A', 



JV f . 



k k± ,^2 ,&4 

d%Xk 1 Xk 2 Xk 3 Xk 4 $k 3 +k 4 ,k 1 +k 2 

. ao hup 



(8) 



with E in units of 2n^f^f- 



rmJ0 " and 



in units of S|ay, 
N c being the number of terms in the linear combina- 
tion. The first part of the GP functional includes the 
kinetic, harmonic and rotating terms of the functional. 
The above expansion turns out to be similar to the one 
used by Abrikosov in his description of the vortex lattice 



superconductors |22]. The Abrikosov lattice in supercon- 
ducting strips was described using a similar expansion 
by one of the authors 3]. Here we make use of our previ- 
ous experience to minimize semianalytically the GP func- 
tional with respect to the coefficients and the wave vec- 
tors fc, imposing the conservation of particles constraint 



J dV(\^\ 2 -n) =0, 



(9) 



where n is the particle density. For an extended analysis 
of the minimization details see Appendix II. 



VORTEX LATTICES IN EXTREMELY 
ELONGATED BEC 

We begin by addressing the case of purely parabolic 
confinement. An interaction of gN — 1.91 is fixed 
throughout. We minimize analytically when feasible (for 
up to N c — 3), while for a higher number of terms we re- 
sort to numerical minimization. Next we add an infinite 
well to the parabolic confinement to simulate a deviation 
from non-harmonicity. If the infinite well is thinner than 
the lateral extension of the BEC, the bosons feel the ad- 
ditional confinement. The existence of the infinite poten- 
tial allows us to study the BEC for rotation frequencies 
higher than the deconfinement frequency in the spirit of 
Ref.^3 for axisymmetric potentials. According to our 
results presented below, the behavior between elongated 
BEC's trapped by only harmonic potentials and BEC's 
confined by stronger confinement potentials is noticeably 
different even for rotation frequencies smaller than the 
critical one. 
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Parabolic confinement 



In this section we show the behavior of an extremely 
elongated BEC as a function of the ratio ^ or the 
anisotropy A (notice that they are related to maintain 
the condition for critical dcconfincment in the y direc- 
tion). In general terms (which we detail below) we find 
two different families or classes of configurations (solid 
and dashed-dotted lines in Fig. ^| . Those denominated 
" even" from now on have a vortex row in the center of the 
BEC and present even parity with respect to an imag- 
inary line along the main BEC axis. In other words, 
they exhibit mirror symmetry, being invariant under the 
transformation than interchanges right and left. Those 
without a vortex row in the middle (denominated "odd" 
from now on) do not present mirror symmetry, but are 
invariant under a combined left-right interchange and a 
translation by the intervortex distance. Both classes of 
solutions are local minima of the multidimensional en- 
ergetic landscape and evolve independently as the con- 
finement varies. In other words, their different symme- 
try acts as a topological characteristic that keeps the 
two minima apart from each other in the energetic land- 
scape, alternating in the ground state. This topological 
characteristic behaves in this sense as the vorticity in ax- 
isymmetric BEC'sQ or superconducting disksQ- Within 
each class of solutions, vortices enter the BEC progres- 
sive and continuously from the edges as the rotation in- 
creases. In general terms one can see that solutions in 
which the vortices are well-defined can be ground states 
while those where incipient vortices appear on the edge 
of the BEC are typically not ground states (see Figs. El 
and 01). The energy and its derivative are always con- 
tinuous for the odd family with some exceptions for the 
even family that we describe below. Derivative discon- 
tinuities for this family sometimes occur in the ground 
state and add up to the ones that occur when the energy 
of the different classes cross each other (see Fig. 0. At 
the crossing points the ground state changes its topologi- 
cal character by increasing the number of vortex rows by 
one. 

In Fig. 0|we present a zoom of the first crossing point. 
For low rotation frequencies, the BEC ground state is 
just formed by the k = component in the expansion and 
does not contain vortex rows (A in Fig. 0). An excited 
state with a very sparse vortex line can also be found at 
those frequencies (C). In the evolution from a BEC with 
no vortex rows to a BEC with one vortex row, the ground 
state first develops a corrugation (B). This corrugation 
is the precursor of the nucleation at the edge of two vor- 
tex rows (see Fig. |3J) • The states with two vortex rows 
are not the ground state until the vortices are fully de- 
veloped (see below) at higher rotation frequencies. The 
appearance of the corrugated state presents a continuous 
derivative (signaled by the long bar in^J, i.e., the sys- 





FIG. 2: Four examples of configurations found for the even 
family at different values of 1.075, 1.029, 1.021 and 1.008. 
Lengths are in oscillator units (o.l.) and gN = 1.91. The state 
in the third panel does not correspond to a ground state. 





FIG. 3: Same as in Fig. 2 for the odd family at different 
values of ^: 1.075, 1.021, 1.008 and 1.006. The first and the 
third states do not correspond to ground states. 



tern presents a second-order phase transition before the 
row of vortices suddenly penetrate when the energy of 
the even family becomes lower and the solution with one 
vortex row becomes the ground state. The second-order 
phase transtion is not surprising since the symmetry that 
characterizes the odd class does not change as the corru- 
gation developes. The corrugated states shown in Fig. [5] 
can be formed by up to five terms. The appropiate varia- 
tional flexibility is here crucial for the correct description 
of this transition. 

As mentioned, derivative discontinuties not only occur 
at the crossing points where the ground state changes 
its number of rows and its symmetry. Once the vortex 
row has appeared in the ground state, the sparse linear 
density of vortices becomes suddenly more dense (D) in 
a first-order phase transition (short bar in Fig. 0J>. As 
far as we have been able to check in our numerics, this 
unexpected sudden change in the ground state is not an 
numerical artifact, although we cannot discard that a 
more complete expansion of the GP wave-function could 
smooth this transition into a second-order one. 

The entrance of a second vortex row at the second 
crossing point also presents some peculiarities. As seen 
in Fig. El before the crossing with the odd family oc- 
curs again, a solution with one vortex row, but differ- 
ent symmetry from the ground state (E), becomes the 
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FIG. 4: Transition from no rows (solid line) to one vortex 
row (dashed-dotted line). Interaction strength is gN = 1.91. 
Additional lines separate different configurations within the 
same family. The corresponding densities are shown in Fig. |S| 



FIG. 6: Same as in Fig.|3] but for the transition from one row 
(dashed-dotted line) to two vortex rows (solid line). Notice 
the crossing intermediate state (dashed line, N c — 5). 
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FIG. 7: BEC densities for the points E, F, and G in Fig. © 
Lengths are in oscillator units (o.L). 



FIG. 5: Panels A, B, C, and D correspond to BEC densities 
at the points shown in|I] Lengths are in oscillator units (o.L). 
The transition from zero (A, N c — 1) to one vortex row (D, 
N c = 2) presents two intermediate type of states (B, N c = 3 
and C, N c = 4). 



ground state (F). This new solution has a snake-like form. 
The fact that the symmetry is different explains why the 
transition is of first order. This solution is, however, al- 
most inmediately overcome by the two-row solution (G). 
From there on the number of rows increases one by one as 
shown in Figs. HG1 and|3 Numerical limitations prevent 
us from making a thourough study of the crossing points 
as we have done in the first two, but we expect similar 
intermediate solutions to appear near the crossing points. 



Parabolic + Infinite well confinement 

Fixing the interaction energy gN = 1.91 and including 
an infinite well of width 7ao we find a similar zero-one 
vortex row transition behavior for a non-parabolic con- 
finement. This is expected since the dispersion relations 
for the parabolic and quasi-parabolic cases are fairly sim- 
ilar in the center of the condensate (for the chosen width) . 
For a small rotation frequency the vectors k of the ex- 
pansion place themselves near the center and are nearly 
the same as those for the parabolic case. However for 
narrower wells, i. e., 4ao, the dispersion relation deviates 
strongly from the parabolic case due to the presence of 
the infinite potential. For the 7ao system we find that the 
one-two vortex row transition does not show the snake- 
like intermediate ground state. This solution has always 
a greater energy than the ground states with two vortex 



FIG. 8: Transition from one to two vortex rows for a quasi- 
parabolic potential defined by a parabolic and an infinite well 
of width 7ao- The zig-zagging configuration is not a ground 
state for any rotation frequency. The correspondent densities 
are similar to the densities in Fig. 




FIG. 9: Transition from zero to one vortex row for a quasi- 
parabolic confimente with an infinite potential of width 4ao- 
All the intermediate states dissapear from the transition. 



rows as shown in Fig. [5] 

For the well of width 4ao the transition from no vor- 
tex to one vortex row losses all the intermediate ground 
states (see Fig. |§J. These results lead us to conclude 
that the infinite potential suppresses intermediate config- 
urations between transitions that changes the number of 
rows by one. We also find that the effects of the boundary 
in the condensate favors configurations without regular 
triangular form. For example, when the condensate rota- 
tion frequency is close to the critical frequency, the extra 
confinement induces states with a quasi-square lattice of 
vortices with energies near the triangular ground state. 
For pure parabolic confinement this state does not appear 
in our calculations. Although a square lattice was found 



FIG. 10: The effects of the infinite confinement produces that 
non triangular configurations appears. However the square 
lattice is not a ground state of the system. 




go 




FIG. 11: Densities for the condensate with an infinite well of 
width 7ao, gN = 1.91, and ^ — 1.001. Lengths are in oscil- 
lator units (o.L). These configurations correspond to those in 
Fig. EH 



experimentally in non-equilibrium elongated svstems|l3j]. 
here, the square lattice only appears when an aditional 
confinement potential is added to the condensate and it 
does not seem to become the ground state of the system 
for as far as we have been able to check. The functional 
energies and the densities of these states are shown in 
Figs. [El and EU 

For rotating frequencies greater than coq the trapping 
potential has the form of an inverted parabola. How- 
ever, for frequencies close to the critical frequency the 
triangular vortex lattice is still the ground state. Finally 
for ultrafast rotation the BEC splits into two channels 
at the edges. This is shown in Fig. ^| for a 7ao infi- 
nite potential and frequencies ^=0.99, 0.98, 0.97 and 
0.96. As shown, this transition presents strong similar- 
ities with the second-order phase transition that takes 
place in superconducting strips at the onset of surface 
superconductivity 23] . 



CONCLUSIONS 

We have minimized the Gross-Pitaevskii functional 
with the usual particle number constraint for the prob- 
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APPENDIX I 

The differential equation for the expansion functions 
Xk is 



dx 2 



Xk 



>2Rx 



2 Q 
1= — ' 



Xk=0, (10) 




FIG. 12: BEC densities for ^=0.99, 0.98, 0.97, and 0.96. 
Lengths are in oscillator units (o.l.) and gN = 1.91. 



lem of a BEC in an infinitely long rotating channel. We 
find the appearance of vortex rows in the BEC along the 
longitudinal direction as the rotation increases. The so- 
lutions alternate between those with a vortex row in the 
center of the channel and those with a maximum in the 
BEC density. For simple parabolic confinement we find 
intermediate states with non-trivial topologies. When an 
infinite potential is added at the boundary of the BEC, 
this intermediate states do not appear in the ground state 
anymore. This extra confinement produces higher energy 
solutions in the form of a square vortex lattice and for 
ultrafast rotation the BEC splits into two channels at the 
edges. 



with et = efe — k 2 



efc the eigenenergies, and 



R = 1 + iV. Making a change of variable u = 
a (\/2Rx S=—k) , we rewrite the equation as 

d 2 



k a 



Xk = 0, 



(11) 



where a is a constant. Choosing a 
reduces to 



4 jr, the equation 



du 2 



Xk 



'2R 



Xk = 0. 



(12) 



The last equation is similar to the quantum harmonic 
oscillator equation, so we identify the analytical solutions 
for the eigenfuctions and the eigenenergies easily: 



and 



Xk = Ce " 2 H n (u) 



e* k = {2n + l)V2R, 



(13) 



(14) 



with H n the Hermite polynomials of grade n and C a 
normalization constant. In the lowest Landau level, (n — 
0), the solutions are finally 
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and 



e k = Vm + k 2 




(16) 



As we expected the dispersion relation shows that if the 
rotating frequency approaches the harmonic frequency 
(Q « ljq) the parabolic confinement dissapears for all the 
wavevectors k. 
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APPENDIX II 



For one component (N c = 1) we have a straightforward 
minimization process. The functional is 



E = 6k\Ck\ 2 + ||Cfc| 4 /fc 



(17) 



with I k — J x\dx. Imposing the number of atoms re- 
striction 

E = {e k -^)\C k \ 2 + 9 -\C k \ A I k +^N (18) 

and minimizing with respect to the coefficient C k we ob- 
tain 



and 



with N 



C k = VN, 



E m in = £kN + -N I k , 



h i = e k + gNI k , 



(19) 



(20) 



(21) 



N„ 



S and N p is the number of particles. For 

this case the ground state corresponds to k = 0. We 
repeat this process for several values of N c and we select 
the lowest E min as ground state. 

For a two-term linear combination (N c = 2) the con- 
straint is TV = 2 1 Cfc | 2 . Using the symmetry we suppose 
that the state which minimizes the functional is com- 
posed by the states (k,—k) and the coefficients C k and 
C_fc are equal. According to these assumptions 

E = 2(e k -fi)C 2 k +g\C k \ 4 A k +^N, (22) 

where A k = I k + 2I k _ fe and I k _ k = J dxxlX-k- Us- 
ing the number of particles constraint and minimizing 
respect the coefficients C k we find 



En 



e k N- 



(23) 



The solutions for the coefficient and the chemical poten- 
tial are 



and 



C k = 



N 
2"' 



lA k N. 



(24) 



(25) 



Finding the value of k that minimizes E m i n , we can es- 
tablish the the density of the BEC. In this case the BEC 
profile is one centered vortex row. The last possible 
analytical calculus is with N c ~ 3. For this configu- 
ration we consider a solution (fc,0, — k) with the same 



value for the coefficients C' k and C'- k as in N c = 2 
case, but including a relative phase for the coefficient 
Co = |Co|e^°. The energy functional with the imposed 
constraint (2C\ + |C | 2 = N) is 

E = 2(6 fe - M )C fe 2 + (6 - Ai )|Co| 2 + | 

{2A k C i k + Jo|Co| 4 + 4S k ^ C 2 \C \ 2 ) + »N(26) 

with S ki< f, = 2I kt Q + cos(2c/)o)Ik.o,-k and where a new 
type of overlap integral between the states, Ik,o,-k = 
J dxxkX-kXoi appears and represents the correlation be- 
tween the vortex rows in the condensate. The phase 
0o which minimizes the functional is 4>o — -|. Calling 
S k = 2I kt0 — I k fi,-k the three-term functional becomes 

E = 2( efc - M )C fe 2 + (£0-^)10) | 2 + § 

(2A k C£ + /(0)|Co| 4 + 4S k C 2 \C \ 2 ) + nN.(27) 

Minimizing the functional with respect the coefficients 
we find 



2 _ e k -e +gN(S k -I(0)) 
g(4S k -2I(0)-A k ) ' 



ct = 



|C | 2 = 



2 _ 2e k - 2e + gN(2S k 



A k ) 



g(4S k - 21(0) - A k ) 



(28) 



(29) 



= 2e k (S k - 1(0)) + e (2S k - A k ) + gN(2S% - 1(0) A k ) 
M 4S k - 2/(0) - A k 

(30) 

and the minimized functional becomes 

(efc -e ) 2 
mm ~ g(AS k ~2I(0)-A k ) + 

N(2e k [S k -I(0)} + e [2S k -A k }) 

4S k - 21(0) - A k 
gN 2 (2Sl-I(0)A k ) 
2(4S k -2I(0)-A k )- 

As far as the chemical potential is concerned 
BE ■ 



+ 



dN 



= A*, 



is easy to see that satifies the expression 

E k in + Eho + E rot + 2Ei nt 



/i = 



N 



(31) 



(32) 



(33) 



found after the direct integration of the GP equation, 
where E k i n , Eh Ql E rot and Ei nt are the kinetic, harmonic 
oscillator, rotational and interaction energies. For three 
components 

E k in + Eho + E rot — 2e fe C 2 + eo|C | 2 , (34) 
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and the interaction energy is 

E int = 9 -{2A k Ci + /(0)|C | 4 + AS k C 2 k \C Q \ 2 ) (35) 

Replacing the analytical coefficients we recover the chem- 
ical potential expression found from the minimization. 
We try minimize the functional for 4 o more terms semi- 
analytically but our results only match those obtained 
numerically when the overlap between terms produces 
vortex rows. We think that this happens because of the 
confinement trapping that maintains the states k close 
to each other. Thus we perform a numerical minimiza- 
tion keeping our restrictions about the behavior of the 
coefficients and the phases that relate them. 
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